• Еще раз про множественные сравнения
• Пример: Поведение песчанок в тесте открытое поле
• perMANOVA
• Условия применимости perMANOVA
• Более подробная интерпретация результатов perMANOVA
• Более сложные дизайны в perMANOVA

Вы сможете

• Вспомнить основные дизайны для тестирования гипотез в рамках дисперсионного анализа
• Применить функции R для тестирования гипотез с помощью perMANOVA

\(\alpha\) — это вероятность совершить ошибку первого рода при тестировании гипотезы (= вероятность отвергнуть истинную нулевую гипотезу, = вероятность найти различия там, где их нет).

Обычно принимается, что \(H_0\) отвергают на уровне значимости \(\alpha = 0.05\).

Когда у нас два средних — все просто, сравнение всего одно. Естественно, вероятность совершить ошибку I рода для группы сравнений \(\alpha _{familywise}\) равна уровню значимости для единственного сравнения \(\alpha _{per\,comparison}\).

\(\alpha _{familywise} = \alpha _{per\,comparison}\)

Но если сравнений много, то растет вероятность совершить хотя бы одну ошибку I рода (найти различия там, где их нет).

Например, если мы хотим попарно сравнить три значения, нам понадобится сделать 3 сравнения.

Пусть мы решили, что в каждом из сравнений будем использовать уровень значимости \(0.05\).

Тогда в каждом из сравнений вероятность совершить ошибку первого рода будет \(\alpha_{per\,comparison} = 0.05\).

Если сделать \(N\) тестов, то вероятность совершить хотя бы одну ошибку I рода в группе тестов (family-wise error rate, FWER) значительно возрастает.

\[\alpha_{familywise} = 1 - (1 - \alpha_{per\,comparison})^N\]

\[\alpha_{familywise} = 1 - (1 - \alpha_{per\,comparison})^N\]

В таблице даны значения \(\alpha _{familywise}\) для разного числа сравнений, если \(\alpha _{per\,comparison} = 0.05\):

1. Взять более жесткий порог уровня значимости. (Например, ввести поправку Бонферрони, поправку Хольма-Бонферрони, воспользоваться процедурой Бенъямини-Хохберга).
2. Изменить схему тестирования гипотезы — тестировать не три независимых гипотезы, а одну сложную (так это, например, происходит в ANOVA).

\(\alpha^*_{per\,comparison} = \alpha_{familywise} /{N}\)

Это жесткий способ, т.к с возрастанием числа сравнений резко снижается уровень значимости и мощность теста. В результате растет риск не найти различий, где они есть.

Ниже \(\alpha _{per\,comparison}\) после поправки Бонферрони, сохраняющие \(\alpha _{familywise} = 0.05\):

Метод Хольма-Бонферрони — это пошаговая процедура.

Чтобы зафиксировать \(FWER \le \alpha_{familywise}\):

1. Сортируем в порядке возрастания \(N\) значений \(p\), полученные в тестах, и присваиваем каждой ранг \(j\) от 1 до \(N\):
   \[p_{1} \le p_{2} \le \cdots \le p_{N - 1} \le p_{N}\]
2. Вводим поправку для значений \(p\)
   \[{p^*_{j}} = min{\{(N - j + 1) \cdot p_{j},\;1\}}\]
3. Сравниваем с \(\alpha\). Отвергаем \(H_0\), если \(p^*_{j} < \alpha_{familywise}\)

В таблице ниже даны результаты нескольких сравнений (\(N = 5\)). С помощью поправки Хольма-Бонферрони для каждого \(p_j\) мы получим свое скорректированное значение \(p^*_{j}\).

\(SS_{total} = SS_{between} + SS_{within}\)

Для этого сравнения Фишер предложил статистику \(F = \frac{SS_{between} / (a - 1)}{SS_{within} / (N - a)}\), где \(a\) — число групп.

F-распределение

Если межгрупповая изменчивость равна внутригрупповой, то \(F\) принадлежит \(F\)-распределению с двумя параметрами \(df_{between} = a - 1\) и \(df_{within} = N - a\), где \(a\) — число классов, \(N\) — общее количество объектов в анализе.

Выявляется влияние фактора А.

Выявляется влияние фактора А, B и их взаимодействия A\(\times\)B.

Некоторые дискретные предикторы могут быть иерархически соподчинены.

Выявляется влияние вложенного фактора B внутри градаций фактора A.

Влияние вложенного фактора B можно оценить, но чаще всего оно не представляет интереса для исследователя.

Влияние блокриующего фактора тоже можно оценить, но часто оно не представляет интереса для исследователя. Обычно блокирующий фактор рассматривается как случайный.

Гипотеза: Разные виды песчанок различаются по поведению в тесте “Открытое поле”

Виды:

• Карликовая песчанка (Gerbillus gerbillus)
• Монгольская песчанка (Meriones unguiculatus)
• Жирнохвостая песчанка (Pachyuromys duprasi)

Фото: (1) XV8-Crisis; (2) Alastair Rae; (3) P.H.J. (Peter) Maas / www.petermaas.nl

Гипотеза: Разные виды песчанок различаются по поведению в тесте “Открытое поле”

Оценка поведения песчанок трех видов по семи признакам:

• Время до выхода в квадрат открытого поля
• Количество актов мочеиспускания
• Количество актов дефекации
• Количество пересеченных квадратов
• Число вертикальных стоек
• Количество актов смещенной активности
• Время проведенное в центре квадрата открытого поля

Фото: (1) XV8-Crisis; (2) Alastair Rae; (3) P.H.J. (Peter) Maas / www.petermaas.nl

pesch <- read.csv("data/pesch.csv", header = TRUE, sep = ";")
head(pesch)

• Какую меру различия можно использовать с этими данными?
• Как можно преобразовать данные?
• Постройте ординацию объектов в осях MDS и раскрасьте точки в соответствии с видами

Поскольку измеренные признаки варьируют в разных масштабах, целесообразно логарифмировать данные.

options(digits = 4)
log_pesch <- pesch
log_pesch[, 3:ncol(pesch)] <- log(pesch[, 3:ncol(pesch)] + 1)
head(log_pesch)

##   Gender Species Time_to_entrance Urination Defecation Quadr_Number
## 1      f    karl            0.000    0.6931      0.000        3.871
## 2      f    karl            3.045    0.0000      1.386        5.762
## 3      m    karl            5.204    0.0000      0.000        5.182
## 4      f    karl            0.000    0.0000      0.000        3.497
## 5      f    karl            4.942    0.0000      0.000        5.328
## 6      m    karl            0.000    0.0000      1.099        3.664
##   Vert_Number Displ_Act Time_in_Centre
## 1       2.485     1.609          0.000
## 2       4.078     1.386          1.946
## 3       4.025     1.946          1.099
## 4       1.792     1.386          0.000
## 5       3.401     2.398          1.386
## 6       2.303     2.197          0.000

library(ggplot2)
library(vegan)
theme_set(theme_bw())

mds_pesch <- metaMDS(log_pesch[, 3:ncol(pesch)], distance = "euclidean")
mds_pesch <- as.data.frame(mds_pesch$points)
mds_pesch$Species <- pesch$Species

ggplot(mds_pesch, aes(x = MDS1, y = MDS2, colour = Species)) + 
  geom_point(size = 5)

Каким способом можно ответить на этот вопрос?

• Мы могли бы взять каждый из поведенческих признаков отдельно и провести одномерный однофакторный дисперсионный анализ.
• Но нас интересует поведение в целом — нужен многомерный анализ.

Что делать если мы хотим оценивать объект по многим признакам сразу?

Примеры:

• Сообщество как целое (много видов)
• Поведение как целое (много элементов)
• Метаболом как целое (много метаболитов)
• Ответы респондентов на множество вопросов в анкетах

Варианты решений:

• MANOVA (Fisher, 1925, Wilks, 1932)
• distance-based Redundancy Analysis (db-RDA) (Legendre, Anderson, 1999)
• perMANOVA (Anderson, 2001; McArdle, Anderson, 2001)

Давно разработан параметрический метод MANOVA (Multivariate Analysis Of Variance). Он дает возможность проводить анализ аналогичный ANOVA. В основе MANOVA лежат представление о многомерном нормальном распределении и расстояниях между центроидами.

В MANOVA сравниваются:

• отклонения точек от групповых центроидов (аналог \(SS_{within}\))
• отклонения групповых центроидов от общего центроида (аналог \(SS_{between}\)).

Сумма квадратов отклонений от центроидов равна сумме квадратов взаимных расстояний, деленной на число объектов

Для Евклидовых расстояний эта закономерность была известна давно (например, Kendall, Stuart 1963).

Anderson, 2001

Anderson, 2001

\(SS_{total} = \frac{1}{N}\sum \sum {d_{ij}^2}\)
Сумма квадратов взаимных расстояний — это сумма квадратов субдиагональных элементов, деленная на число объектов N.

\(SS_{within} = \frac{1}{n}\sum \sum {d_{ij}^2 \cdot \epsilon_{ij}}\)
Внутригрупповая сумма квадратов — это сумма всех сумм квадратов расстояний между элементами для каждой группы, деленная на n число объектов в группе.

Тогда межгрупповая сумма квадратов \(SS_{between} = SS_{total} - SS_{within}\)

\[F = \frac{SS_{between} / (a - 1)}{SS_{within}/(N - a)}\]

Для оценки значимости псевдо-F используется пермутационная процедура:

• Случайным образом перетасовываются строки исходной матрицы
• После каждого акта пермутации вычисляется \(F_{perm}\)
• Уровень значимости
• Внимание! Для иерархического дизайна процедура пермутации имеет свои особенности (обсудим позднее).

\[p = \frac {N_{F_{perm} \ge F}} {N_{permutations}}\]

library(vegan)
permanova_pesch <- adonis(log_pesch[3:9] ~ Species, data = log_pesch,
                          method = "euclidean")
permanova_pesch

## 
## Call:
## adonis(formula = log_pesch[3:9] ~ Species, data = log_pesch,      method = "euclidean") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)  
## Species    2        60   30.02    3.58 0.237  0.012 *
## Residuals 23       193    8.39         0.763         
## Total     25       253                 1.000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Мы видим традиционную для ANOVA таблицу результатов. Что здесь что?

## 
## Call:
## adonis(formula = log_pesch[3:9] ~ Species, data = log_pesch,      method = "euclidean") 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model    R2 Pr(>F)  
## Species    2        60   30.02    3.58 0.237  0.012 *
## Residuals 23       193    8.39         0.763         
## Total     25       253                 1.000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Поведение разных видов песчанок в тесте открытое поле значимо различалось (perMANOVA, p = 0.012).

1. Требуется равенство внутригрупповых дисперсий (гомоскедастичность).

2. Желательно использование сбалансированных данных (с равным объемом групп), т.к. этом случае perMANOVA устойчив к гетерогенности дисперсии (Anderson, Walsh, 2013).

Для проверки можно использовать функцию betadisper(), изначально предназначенную для сравнения \(\beta\)-разнообразия сообществ.

Функция betadisper() вычисляет внутригрупповые центроиды и координаты точек в пространстве главных координат (Principal Coordinate Analysis = PCoA = metric MDS).

Затем, значимость различий отклонений от центроидов в разных группах проверяется с помощью процедуры PERMDISP2.

dist_pesch <- vegdist(log_pesch[,3:ncol(pesch)], method  = "euclidean")
PCO_pesch <- betadisper(dist_pesch, log_pesch$Species)

Объект, возвращаемый betadisper(), позволяет также нарисовать наши объекты в пространстве главных координат (PCoA).

plot(PCO_pesch, main = "PCoA ordination")

Процедура PERMDISP2 реализована в пакете vegan в функции anova.betadisper(). Это многомерный аналог теста Левина на гомогенность дисперсий в группах, который иногда используется для проверки условий применимости дисперсионного анализа.

anova(PCO_pesch)

## Analysis of Variance Table
## 
## Response: Distances
##           Df Sum Sq Mean Sq F value Pr(>F)
## Groups     2    9.7    4.86    2.05   0.15
## Residuals 23   54.5    2.37

• Не выявлено значимых различий разброса внутригрупповых расстояний.

boxplot(PCO_pesch)

На приведенной ординации видно, что точки, соответствующие Монгольским песчанкам расположены отдельно от остальных.

Для выявления попарных различий нужны попарные сравнения.

Внимание! В пакете vegan пост хок тест не реализован. Но мы можем сделать простейшую его версию самостоятельно.

Проведем попарные сравнения между группами, то есть:

• Карликовые vs. Монгольские
• Карликовые vs. Жирнохвостые
• Монгольские vs. Жирнохвостые

Фото: (1) XV8-Crisis; (2) Alastair Rae; (3) P.H.J. (Peter) Maas / www.petermaas.nl

pairwise_permanova <- function(dat, group, strata = NULL, ...){
  pair <- combn(unique(as.character(group)), 2)
  ncomb <- ncol(pair)
  res <- rep(NA, ncomb)
  for (i in 1:ncomb) {
    filter <- group %in% pair[, i]
    if(is.null(strata)){
      posthoc <- adonis(dat[filter, ] ~ group[filter], ...)$aov.tab$Pr[1]
    } else {
      posthoc <- adonis(dat[filter, ] ~ group[filter], 
                        strata = strata[filter], ...)$aov.tab$Pr[1]
    }
    res[i] <- posthoc
    names(res)[i] <- paste(pair[, i], collapse = " vs. ")
  }
  return(res)
}

p_vals <- pairwise_permanova(
  dat = log_pesch[, -c(1:2)], group = log_pesch$Species, 
  method = "euclidean", permutations=9999)
p_vals

##         karl vs. mongol   karl vs. zhirnokhvost mongol vs. zhirnokhvost 
##                  0.0011                  0.4029                  0.0031

Это все? Пишем статью?

Фото: (1) XV8-Crisis; (2) Alastair Rae; (3) P.H.J. (Peter) Maas / www.petermaas.nl

Мы делали три пары сравнений — нужно ввести поправку для уровня значимости. Будем считать значимыми различия в тех сравнениях, где после введения поправки \(p < 0.05\).

Можно посчитать скорректированные значения уровня значимости \(p\) с учетом поправки Хольма-Бонферрони.

p.adjust(p_vals, method = "holm")

##         karl vs. mongol   karl vs. zhirnokhvost mongol vs. zhirnokhvost 
##                  0.0033                  0.4029                  0.0062

Выясним, влияет ли пол и вид песчанок на поведение.

Отфильтруем исходные данные (в случае с жирнохвостыми песчанками были изучены только самки)

log_pesch2 <- log_pesch[log_pesch$Species != "zhirnokhvost", ]

• Создайте датафрейм из файла simulated_data.csv (Это данные симулированные по алгоритму, приведенному в справке по функции adonis())
• В этом датафрейме записано обилие двух видов на экспериментальных площадках двух типов: без добавления и с добавлением NO3, по 6 повторностей в каждом эксперименте.
  Эксперименты были независимо проведены на 3 полях.
• Оцените, зависит ли структура совместного поселения этих двух видов от концентрации NO3.

com <- read.csv("data/simulated_data.csv", sep = ',', header = T)

# Ошибочный дизай
com_permanova <- adonis(com[,1:2] ~ com$NO3)

# Правильный дизайн
com_permanova2 <- adonis(com[,1:2] ~ com$NO3, strata = com$field)

• При одновременном тестировании нескольких гипотез растет вероятность ошибки I рода.
• Чтобы контролировать вероятность ошибки I рода, либо используют более жесткий уровень значимости, либо тестируют сложную гипотезу вместо нескольких простых.
• perMANOVA дает возможность тестировать сложные гипотезы в отношении явлений, описанных по многим переменным (т.е. на многомерных данных).
• В perMANOVA можно использовать любые коэффициенты различия.
• Для применения perMANOVA требуется равенство разбросов точек между центроидами их групп, но при равных объемах групп анализ устойчив к отклонениям от этого условия.
• При использовании perMANOVA важно не запутаться в дизайне.

• Primer 6.0 + PERMANOVA Коммерческий продукт.
• PAST Некоммерческая программа. Здесь метод называется NPMANOVA.
• Оригинальная программа М. Андерсон PERMANOVA и PERMDISP.

• Anderson, M.J. 2001. A new method for non-parametric multivariate analysis of variance. Austral Ecology, 26: 32–46.
• Anderson, M.J. 2005. PERMANOVA: a FORTRAN computer program for permutational multivariate analysis of variance. Department of Statistics, University of Auckland, New Zealand.
• Anderson, M.J. (2004). PERMDISP: a FORTRAN computer program for permutational analysis of multivariate dispersions (for any two-factor ANOVA design) using permutation tests. Department of Statistics, University of Auckland, New Zealand.
• Legendre P., Legendre L. (2012) Numerical ecology. Second english edition. Elsevier, Amsterdam.